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FOREWORD 


This  report  describes  work  done  by  Sabbagh  Asso^ates,  Inc.  under  contract  N60921- 
85-C-0046  with  the  Naval  Surface  Weapons  Center,  White  Oak,  as  part  of  the  Center’s 
Small  Business  Innovation  Research  Program  (SBIR)»The  research  problem  is  to  develop 
a  model  and  an  inversion  algorithm ‘^haH^  suitable  for  the  three-dimensional  quantitative 
nondestructive  evaluation  (NDE)  of  advanced  composite  materials  by  using  eddyicurrents. 
The^cfcpproach  is  based  on^^abbagh  Associates’  work  in  eddy-current  NDE  of  conventional 
metals.  The^echnical  objectives  are  to  determine  the  feasibility  of  using  multifrequencies 
for  this  job,  to  determine  in  localized  regions  the  fiber-resin  ratio  in  graphite  epoxy,  and 
to  determine  more  precisely  the  types  of^omalies,  whether  flaws,  delaminations,  broken 
fibers,  etc.,  that  can  be  reconstructed  b^tve^r  inversion  method.  These  objectives  are  met 
by:  (1)  applying  rigorous  electromagnetic  theory  to  determine  a  Green’s  function  for  a 
slab  of  anisotropic  composite  material,  (2)  determining  the  integral  relations  for  the  direct 
and  inverse  problem^^  using  the  Green’s  function  just  derived,  (3)  determining  suitable 
numerical  algorithms  for  solving  the  inverse  problem,  and  (4)  writing  a  computer  program 
to  execute  the  model.  These  objectives  have  been  met;  in  addition  we  illustrate  some 
interesting  electromagnetic  phenomena  in  composite  materials.  *  !■ 


Composite  materials  in  the  form  of  fiber-reinforced  matrix  materials  ais,  for  example, 
graphite-epoxy,  are  being  increasingly  used  in  critical  structures  and  structural  compo¬ 
nents  because  of  their  high  strength-to-weight  ratio.,  Such  structures  range  from  high- 
performance  aircraift  to  rocket  motor  cases  to  consumer  goods,  such  as  automobiles,  golf 
clubs  and  tennis  rackets.  In  order  to  assess  the  integrity  of  these  structures,  thereby  de¬ 
termining  the  presence  and  size  of  flaws,  it  is  necessary  to  employ  suitable  methods  of 
quantitative  nondestructive  evaluation,  such  2is  eddy-currents.  This  research  will  have 
significant  applications  in  the  military  and  commercial  aircraft  industries,  as  well  as  auto¬ 
mobiles  and  consumer  products. 


Approved  by: 


JACK  R.  DIXON,  Head 
Materials  Branch 
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CHAPTER  1 

IDENTIFICATION  AND  SIGNIFICANCE  OF  THE  PROBLEM 

The  proposed  research  problem  is  to  develop  a  model  and  an  inversion  algorithm  that 
is  suitable  for  the  quantitative  nondestructive  evaluation  (NDE)  of  advanced  composite 
materials  by  using  electromagnetic  methods,  especially  eddy-currents.  We  are  specifically 
interested  in  determining,  in  localized  regions,  the  fiber-resin  ratio  in  graphite-epoxy. 

Composite  materials  in  the  form  of  fiber-reinforced  matrix  materials  as,  for  example, 
graphite-epoxy,  are  being  increasingly  used  in  critical  structures  and  structural  components 
because  of  their  high  strength-to-weight  ratio.  In  order  to  assess  the  integrity  of  these  struc¬ 
tures,  it  is  necessary  to  employ  suitable  methods  for  quantitative  NDE.  One  method  uses 
eddy-currents;  composite  materials,  however,  are  inherently  anisotropic,  which  means  that 
many  of  the  classical  eddy-current  technology  and  design  procedures  are  not  applicable.  In 
addition,  composite  materials  vary  widely  in  their  permittivities  zmd  conductivities,  which 
means  that  new  analyses  must  be  carried  out  to  develop  effective  strategies  for  using  eddy- 
currents  in  quantitative  NDE.  A  final  problem  is  that  there  is  a  variety  of  potential  failure 
modes  in  composites,  such  as  delaminations,  fiber-breakage  due  to  impact  damage,  flaws, 
etc.,  some  of  which  may  not  be  readily  detectable  by  eddy-currents  [l].  In  order  to  com¬ 
plement  the  empirical  studies  of  [1]  it  is  necessary  to  embark  upon  a  rigorous  quantitative 
NDE  program  for  composites  to  assess  the  role  that  eddy-currents  play  in  it,  and  especially 
to  determine  a  suitable  inversion  algorithm. 

In  quantitative  NDE  one  seeks  to  determine  the  orientation  and  spatial  extent  of  flaws, 

i.e.,  to  reconstruct  or  image  flaws.  By  so  quantifying  a  flaw  one  is  better  able  to  assess  its 
potential  for  damaging  the  structure  in  which  it  is  located,  thereby  offering  a  means  for 
judging  whether  to  eliminate  the  flawed  member  or  not.  Clearly,  such  decisions  can  save 
hundreds  of  thousands  of  dollars  in  “down  time”,  as  well  as  potentially  saving  lives. 

The  present  state  of  the  art  in  quantitative  eddy-current  NDE  of  composites  does  not 
include  the  ability  to  reconstruct  or  image  anomalous  regions,  or  flaws,  in  three  dimensions. 
The  objective  of  Phase  I  research,  therefore,  is  to  develop  a  model  and  algorithm  suitable 
for  inverting  eddy-current  data  for  the  reconstruction  of  flaws  in  quantitative  NDE  of 
advanced  composite  materials.  The  specific  technical  objectives  are: 

1.  To  develop  a  model  describing  the  interaction  of  induced  eddy-currents  and  the 
conducting  graphite  fibers  in  graphite  epoxy, 

2.  To  develop  a  fiber  density  algorithm  based  on  this  model, 

3.  To  determine  the  feasibility  of  using  multifrequencies  for  this  job. 

These  objectives  can  be  met  by  accomplishing  the  following: 

1.  Apply  rigorous  electromagnetic  theory  to  determine  a  Green’s  function  for  a 
stratified  half-space,  or  a  finite  slab,  of  anisotropic  composite  materials. 
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2.  Determine  the  integral  relations  for  the  direct  and  inverse  problem,  using  the 
Green’s  function  derived  in  1. 

3.  Determine  suitable  numerical  algorithms  for  solving  the  inverse  problem. 

4.  Write  a  computer  program  to  execute  the  model. 

Electromagnetic  theory  and  the  Green’s  function  are  dealt  with  in  Chapters  2-4, 
and  in  Appendices  A,  B  and  C.  Integral  equations  for  the  direct  and  inverse  problems 
are  derived  in  Chapter  5,  which  also  includes  a  discussion  of  the  discretization  of  these 
equations  by  means  of  the  method  of  moments  and  the  use  of  multifrequencies.  Chapter 
6  describes  the  numerical  algorithms  that  were  used:  these  are  least-squares,  by  means 
of  the  singular  value  decomposition  and  QR-decomposition,  together  with  the  Levenberg- 
Marquardt  stabilization  parameter,  and  an  iterative  algorithm  for  analytic  continuation 
(or  Fourier  extrapolation.) 

We  conclude  the  report  with  a  description  of  some  numerical  experiments,  comments 
and  conclusions. 
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CHAPTER  2 

AN  ELECTROMAGNETIC  MODEL  FOR  GRAPHITE-EPOXY 
BACKGROUND 

Eddy-current  methods  for  the  examination  of  carbon  fiber  reinforced  epoxy  resins 
and  other  composite  materials  have  been  discussed  and  analyzed  by  Owston  and  Prakeish 
[2]-[5].  These  analyses  have  been  based  on  an  ad  hoc  equivalent  circuit  in  which  the 
composite  test  piece  is  regarded  as  being  inductively  coupled  to  the  probe,  much  as  in 
the  classical  treatment  of  eddy-current  evaluation  of  metals.  Though  the  technique  gives 
a  useful  indication  of  the  form  of  the  results,  a  more  satisfactory  approach,  as  Owston  [4] 
points  out,  is  to  use  a  field-theoretic  analysis  which  is  capable  of  giving  exact  results  for  a 
given  model.  A  field-theoretic  analysis  is  also  desirable  when  computing  electromagnetic 
interactions  for  shielding  effectiveness  of  advanced  composites  in  aircraft  [6]. 

The  heart  of  the  problem  is  to  determine  a  Green’s  function  for  the  composite  material. 
(The  Green’s  function  is  the  electromagnetic  field  produced  by  a  point-source  of  current.) 
Much  work  has  been  done  in  recent  years  on  the  subject  of  electromagnetic  interactions 
with  composite  materials,  mostly  in  the  context  of  electromagnetic  shielding  of  avionics 
equipment  from  electromagnetic  pulses  [6j-[10].  Some  of  this  work  is  directly  applicable 
to  the  problem  of  computing  eddy-current  flow  within  composites,  but  the  Green’s  func¬ 
tion  problem  must  be  attacked  by  applying  rigorous  electromagnetic  theory  to  anisotropic 
media. 

CONSTITUTIVE  RELATIONS  FOR  ADVANCED  COMPOSITES 

Advanced  composite  materials  are  laminates  made  up  of  a  number  of  individual  layers 
bonded  together.  Each  layer  consists  of  a  unidirectional  array  of  long  fibers  embedded  in, 
and  firmly  bonded  to,  a  matrix.  The  basic  building  blocks  of  any  specific  composite  are 
defined  by  the  types  of  fibers  and  matrix  involved.  Some  fiber-matrix  systems  are:  boron- 
epoxy,  graphite-epoxy,  Kevlar-epoxy,  graphite-polymide  and  graphite-thermoplaustic  [7]. 

The  matrix  for  each  of  these  materials  is  normally  a  good  dielectric,  whereas  the  fibers 
v^ary  in  electrical  conductivity  from  modest  (graphite)  to  a  poor  dielectric  (boron)  to  a  good 
dielectric  (Kevlar).  These  materials  are  nonmagnetic,  so  that  the  magnetic  permeability 
is  fio* 

Composites  have  anisotropic  conductivities  because  of  the  unidirectional  arrays  of 
fibers  within.  For  example,  for  graphite-epoxy  the  average  macroscopic  conductivity  along 
the  fiber  direction  is  20,000  mhos/m,  whereas  in  the  direction  transverse  to  the  fibers,  the 
conductivity  is  100  mhos/m.  It  may  be  surprising  to  find  a  nonzero  transverse  conductivity 

*  These,  and  other  facts  about  electromagnetic  constitutive  relations,  are  taken  from  [7), 
which  is  the  most  complete  reference  on  electromagnetic  modeling  of  composite  materials 
that  we  have  found. 
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in  graphite-epoxy,  in  view  of  the  earlier  statement  that  the  matrix  is  a  good  dielectric. 
The  fact  is  that  there  is  enough  local  fiber-to-fiber  contact  that  the  average  macroscopic 
conductivity  is  not  zero.  (See  Figure  1(a).)  Other  materials,  of  course,  have  different 
longitudinal  and  transverse  conductivities,  as  shown  in  Table  1  [7]: 

TABLE  1.  SUMMARY  OF  ELECTRICAL  PROPERTIES 
OF  SOME  COMPOSITES 


Graohite-Enoxv 

Boron-Eooxv 

Kevlar 

Permeability  /xr 

1 

1 

1 

Permittivity  £r 

indeterminate 

5.6 

3.6 

DC  Conductivity  (mhos/m) 
longitudinal  cl 

2  X  10“* 

30 

6  X  10-- 

transverse  ot 

100 

2  X  10“® 

6  X  10-- 

Anisotropy  Ratio 

200 

1.5  X  10® 

1 

The  reason  that  cr  for  graphite-epoxy  is  indeterminate  is  because  the  fiber-to-fiber  contact 
effectively  shunts  the  capacitance  between  fibers  with  a  fairly  low  resistance  path,  making 
it  impossible  to  measure  dielectric  permittivities  at  frequencies  less  than  100  megahertz,  or 
so.  Thus,  in  Figure  1(b),  which  shows  a  possible  ac  equivalent  circuit  for  eddy-current  flow, 
the  capacitors  are  effectively  shorted  by  the  fiber-to-fiber  resistors  at  the  lower  frequencies. 

The  anisotropy  of  the  composite  manifests  itself  in  a  complex  permittivity  tensor,  the 
tensor  being  diagonal  in  a  coordinate  system  (^1,^2.  ^s).  where  is  parallel  to  the  average 
fiber  direction,  ^2  is  perpendicular  to  the  average  fiber  direction,  but  lies  in  the  plane  of 
the  composite  layer,  and  $3  is  perpendicular  to  both  fibers  and  the  plane  of  the  layer: 

(en  0  0\ 

I  =  0  C22  0  .  (1) 

\  0  0  633  y 

Here,  £,,  =  with  j  =  V^,  and  u)  is  the  angular  frequency.  The  coordinate 

system  just  defined,  for  which  the  complex  permittivity  tensor  is  diagonal,  is  not  necessarily 
the  laboratory  coordinate  system,  (x,y,z),  in  which  the  electromagnetic  field  vectors  are 
defined.  In  any  case,  the  tensor  symbol  will  be  used,  and  the  components  in  a  particular 
coordinate  system  may  be  computed  by  applying  the  usual  rules  for  transforming  Cartesian 
tensors. 

From  here  on  we  will  consider  only  graphite-epoxy,  for  which  cn  =  £22  =  C33  =  to, 
Oil  =  2  X  10'*  mhos/m,  and  o’22  =  <^33  =  100  mhos/m. 

MATRIX  WAVE  EQUATIONS  FOR  GRAPHITE-EPOXY 

MjLXwell’s  equations  are  the  fundamental  equations  for  electromagnetic  analysis.  Us¬ 
ing  the  constitutive  relations  just  defined  for  graphite-epoxy,  these  equations  become,  in 


the  sinusoidal  steady-state: 
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V  y.  E  —  —juJUoH 

V  X  H  =  j(jjt  •  E  -f 


Because  of  the  anisotropy  of  graphite-epoxy,  it  is  convenient  to  work  with  a  matrix  for¬ 
mulation  of  these  equations  that  has  been  useful  in  crystal  optics,  plasmas  and  microwave 
devices  [ll]-[19].  We  start  by  writing  (2)  in  the  6x6  matrix  form 


(  jwl  -Vx/\  ^ 

\VxI  iunol)\^)  \  0  )' 


where  the  permittivity  tensor  includes  a  conductivity  term,  as  in  (l),  except  that  the 
matrix  is  not  necessarily  diagonal.  I  is  the  identity  matrix,  and  0  is  the  zero  vector. 

We  write  the  matrix  of  partial  derivatives  as: 

(  V  X  7  ”  V  ^ 


where 


U.  = 


0  0  0 

0  0  1 

0-10 


0  0  0 

0  0-1 
0  1  0 


1/2  = 


1/3  = 


0  0  1 

0  0  0 

-10  0 


0-10 
1  0  0 

0  0  0 


0  0-1 
0  0  0 
10  0 


0  10 

-10  0 
0  0  0 


When  these  matrices  are  substituted  into  (3),  along  with  the  definitions 

.  ..  /-/(•)  \  =  /=  s  \ 


(!)■  =  ^=(i  i?)- 


we  get 


juK  -h  Ui  dldx  -f  1/2  djdy  -H  C/3  djdz  e  = 

Upon  defining  the  two-dimensional  Fourier  transforms  in  (x,y)  by 

i{k^,ky,z)  =  ^  jj  e{x,y,z)e^^'^^^+''»^Uxdy 

e{x,y,2:)  =  jj  i{k^,ky,z)e~^^'^^^+’^’>^Uk^  dky, 


(8) (a) 

im 
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we  take  the  Fourier  transform  of  (7); 

^jujK  -  jkxU I  -jkyU2  +  U3d/dz^e  =  -/  \  (9) 

To  get  this  result  we  assume  that  the  composite  material  is  homogeneous  in  the  (x,y) 
plane,  so  that  K  is  independent  of  (x,y). 

The  third  and  sixth  rows  of  (9)  are  independent  of  s-derivatives.  Thus,  we  can  use 
these  two  equations  to  solve  for  C3  and  ce  in  terms  of  the  remaining  variables: 


When  we  use  these  two  equations  to  eliminate  the  z-components  from  the  remaining  equa¬ 
tions  in  (9)  we  end  up  with  system  of  four  first-order  differential  equations  in  the  four 
transverse  electromagnetic  field  components.  This  system  is  written  as  the  4x4  vector- 
matrix  equation: 


det  ,jf 


where 


fit  [41,62,44,65]  —  ^yt  Hxt 
the  subscript,  ‘f’,  denotes  transverse  field  variables. 


The  components  of  S  and  U  are: 


5ii  = 


jkxKzi 


S\2  = 


jkxKz2 


5,3  = 


r 

T7  > 


5i4  =  — iwMO  +  J  — 


531  =  juKzi  —  jw 


523  = 

tjjKzz 

532  =  ~  — 


541  =  —ju}K\\  -t-  jw 


JuKi2 
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CHAPTER  3 

DERIVATION  OF  A  GREEN’S  FUNCTION  FOR  A  GRAPHITE-EPOXY  SLAB 


The  formal  solution  of  (11)  defines  the  Green’s  function,  G: 


lt{z)  =  j  G{zW)  ‘3^  ^{z')dz'. 


The  double  bars  over  G  denote  a  tensor,  and  the  tilde  over  any  variable  means  a  transverse 
Fourier  transform,  as  before.  Clearly,  the  Green’s  function  is  the  transverse  field  at  z  due 
to  a  point  current  source  at  z\  Thus,  it  is  a  4  x  3  tensor,  because  there  are  four  transverse 
field  components,  as  shown  in  (12),  and  the  point  current  source  can  be  oriented  in  the 
y->  oi  z-direction.  Because  the  current  flow  in  the  models  studied  in  this  report  are 
confined  to  the  transverse  plane,  i.e.,  the  (x,y)-plane,  we  will  be  interested  in  only  the  first 

two  columns  of  G. 

The  calculation  of  the  Green’s  function  requires  a  knowledge  of  the  eigenmodes  of  (11) 

for  a  nonstratified  medium.  These  are  solutions  of  the  homogeneous  form  of  (11),  with  5 
constant  with  respect  to  z.  The  eigenmode  theory  is  developed  in  Appendix  A. 

In  order  to  compute  G  we  must  know  the  geometry  of  the  composite  material.  In  this 
report  we  will  work  with  the  plane-parallel  slab  that  is  shown  in  Figure  2,  with  the  x-axis 
parallel  to  the  fiber  direction,  and  the  y-direction  trsmsverse.  Thus,  the  region  of  interest 
consists  of  three  parts:  the  region  above  the  slab  (which  is  free  space),  the  composite  slab, 
and  the  region  below  the  slab.  Call  these  regions  1,  2,  3,  respectively,  as  shown  in  Figure 
2,  and  introduce  the  following  notation  for  the  Green’s  function: 

Gij{z\z')  =Field  produced  at  z  in  region  i,  due  to  a 
point  source  of  current  at  z'  in  region  j ; 
i,j  =  1,2,3. 

As  we  will  see  later,  the  only  Green’s  functions  that  are  needed  are  G21  and  G12. 

In  order  to  compute  G21,  consider  Figure  3(a),  in  which  the  point-source  of  current  is 
located  at  z*  in  region  1.  We  further  subdivide  the  three  regions  into  four  regions,  called 
A-D,  and  expand  the  vector  field  solutions  in  each  of  these  four  regions  in  terms  of  the 
eigenvectors  of  Appendix  A: 


-aio\ 

■«20  i  ft 


-^10  -A„  (*-*') 

1 

—P20  y 


(16) (o) 
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aio\  /— aio\  [  ^  \ 

020  Iao (»-»')  -020  1-Ao (*-*')  4.  e  Lao(*-*') 

I  J  \  1  J  V-/?20/ 

/  0  \ 


+  /[  ^10  e-Ao  (*-*') 

y-p2oJ 


M  i,Ai2_L^/f  <*2  lg-Ai*_j_g/|  1  .Aa*  _L  /•/ (  1  ^-As* 


e^»*  +  /' 


(16)(c) 


gAo(*-l-*o)  ^ 


0 

/Jio 

1 

~y^20 


Ao(ar-t-*o) 


(16) (d) 


where  ai  =  S14/X1,  02  =  524/Ai,  /3i  =  X3/S32,  and  /?2  =  531/532.  The  subscript  ‘0’ 
denotes  quantities  defined  in  isotropic  free-space. 

Equations  16(a)  and  16(d)  contain  only  outgoing  waves  at  ±infinity.  In  order  to 
determine  the  unknown  expansion  coefficients  a  through  h,  c'  through  f\  we  must  satisfy 
certain  boundary  conditions.  The  fields  must  be  continuous  at  the  two  boundaries,  z  = 
0,andz  =  — zo>  which  do  not  contain  current  singularities.  At  z  =  — z',  which  does  contain 
a  point  singularity,  there  must  be  a  discontinuity.  The  amount  of  discontinuity  can  be 
inferred  from  (11),  in  which  the  imposed  source  current  has  the  form: 

A 

where  j  is  the  unit  vector  in  the  direction  of  current  flow. 

Upon  integrating  (11)  an  infinitesimal  distance  across  the  plane  z  =  z',  we  get 

=  U]  (18) 

where  (-f )  denotes  the  limit  immediately  above  z',  and  (— )  denotes  the  limit  immediately 
below. 

In  order  to  find  the  response  to  x-directed  currents  (the  first  column  of  the  Green’s 

A 

tensor)  we  let  j  =  3*.  The  right-hand  side  of  (18)  becomes 


(17) 
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Similarly,  for  a  y-directed  current,  which  will  produce  the  second  column  of  the  Green’s 
tensor,  we  let  j  =  ay,  and  get 

/0\ 


for  the  right-hand  side  of  (18). 


Once  we  obtain  the  boundary  conditions,  the  solution  for  the  expansion  coefficients 
becomes  a  straight-forward  problem  in  linear  algebra. 


The  computation  of  G12  proceeds  in  a  similar  manner,  except  that  the  four  subregions 
are  shifted  downward,  as  shown  in  Figure  3(b).  The  field  expansions  are  given  by: 


— OlO 


A\  a  I 


-Ao  1 
—^20  / 


(21) (a) 


C  “2  + 


-“2  .-A,  (.+•').,  h  ,A,(«+,') 

0  1 


(  0 

+  /  e->3  (*+*')  (21)  (6) 

\-02J 


«2 


~“2  -A, (*-)-*')  ,  Z  P\  -AaC*-!-*') 

0  -re  j  e 


I  °  "1 

+  f'\  f'  (*+*')  (21)  (c) 

\-02J 


J)  :  g  *^20  gAo(*-t-*o)  _j_  ^  ^10  gAo{*-r»o)^ 


—020 


(2l)(d) 


and  the  current  singularity  is  at  z',  which  is  within  the  slab. 


>. 
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The  boundary  conditions,  including  the  discontinuity  condition,  (18),  are  the  same 
as  before.  In  particular,  the  inhomogeneous  terms,  (19)  and  (20),  corresponding  to  an 
x-directed  and  y-directed  current  source,  respectively,  continue  to  hold. 

In  Appendix  B  we  ske^h  the  application  of  the  boundary  conditions  for  solving  for  the 

expansion  coefficients  for  G12.  The  process  for  G2X  is  similar.  After  the  coefficients  have 
been  determined,  we  substitute  c'  through  f  into  (16)  (c)  to  complete  the  computation  of 

G21;  for  Gi2  we  substitute  (a,  6)  into  (21)(a).  These  steps  complete  the  solution  of  the 
Green’s  tensor. 


10 
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CHAPTER  4 

CALCULATION  OF  THE  ELECTROMAGNETIC  FIELD 
WITHIN  A  GRAPHITE-EPOXY  SLAB 

Having  computed  a  Green’s  function  for  the  graphite-epoxy  slab,  it  is  possible  to 

determine  the  field  within  the  slab  due  to  currents  above  the  slab  (by  using  G21),  or  to 

determine  the  field  above  the  slab  due  to  currents  within  the  slab  (by  using  G12).  The 
latter  operation  will  be  used  in  the  inverse  problem.  The  method  of  determining  the  fields, 
given  the  currents,  is  to  use  (15). 

=(*■) 

We  have  considered  four  different  impressed  current  sources,  j  :  an  infinite  current 
sheet,  a  filamentary  circular  current  loop,  and  two  finite  circular  current  sheets  with  dis¬ 
tributed  windings,  one  with  a  uniform  distribution,  and  the  other  with  a  linearly  increasing 
distribution.  Each  of  the  current  sources  is  parallel  to  the  surface  of  the  slab.  The  three 
circular  currents  simulate,  to  a  degree,  flat  ‘pancake’  coils. 

The  infinite  current  sheet  plays  a  significant  role  in  the  inversion  process,  which  will 
be  presented  in  the  next  two  chapters,  and  will  be  discussed  first. 

INFINITE  CURRENT  SHEET 

The  problem  is  illustrated  in  Figure  4.  In  the  Fourier  transform  domain,  the  current 
sheet  is  given  by 

\kx,ky,z')  =  Io6{kx)6{ky)6{z'  —  «")(o*co8  0  -I-  Oy  sin^i),  (22) 

where  Iq  is  the  surface  current  density  in  zunperes  per  meter,  and  $  is  the  direction  of 
current  flow.  The  delta  function  at  the  origin  of  {kx,ky)  space  is  due  to  the  fact  that 
the  current  sheet  is  uniform  over  the  entire  (x,y)  plane.  The  principal  axes  are  the  x- 
and  y-axes,  the  x-axis  being  along  the  fiber  direction,  and  the  y-axis  tremsverse.  We  will 
compute  the  field  within  the  slab  when  the  current  is  oriented  along  either  of  these  axes 
(i.e.,  0  =  0°  or  90°)  ,  and  then  use  a  simple  tensor  transformation  to  compute  the  fields 
when  the  current  sheet  is  oriented  at  any  other  angle. 

=(»■) 

Because  the  spectrum  of  j  is  concentrated  at  the  origin,  we  can  simply  put  kx  — 

Aiy  =  0  in  computing  the  Green’s  function,  Gji.  Of  course,  the  evaluation  of  the  integral 
in  (15)  is  simplified  because  of  the  delta  function  at  z".  These  facts  allow  us  to  write  down 
a  simple  analytic  expression  for  the  transverse  fields  within  the  slab: 


current  in  x-direction: 


(23)  (a) 


m 
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current  in  y-direction: 


where 


Ey 

\HyJ 


gAs*  ^  ji 


.-Ao*" 


9  = 


2aii 


-Iq 


^  (^11 +^o)^e^‘*° -(*?ii ->?o)^e 

{riii+VoV -{vn-rjoV 

_--Ao*" 

/i=— - 7o 

20!22 

_  (»/  +  -  (*?  - 

(»?  +  »7o)^  -  (»?  -  »?o)^ 


Ao  =  jw(/ioCo)^''^,Ai  =  yu>(MoiCn)^^*,A3  =yw(Mo/iC)^''*, 
no  =  (Mo/€o)^^*,nii  =  (Mo/ifn)'/*,f7  =  [lio/Ky^^. 


(23)(b) 


(24) 


We  see  from  (23)  that  an  x-directed  uniform  current  sheet  produces  an  x-directed 
electric  field  and  a  y-directed  magnetic  field,  whereas  a  y-directed  current  sheet  produces 
a  y-directed  electric  field  and  an  x-directed  magnetic  field.  This  is  what  we  expect  from 
an  isotropic  medium;  it  occurs  in  an  anisotropic  medium  when  the  x-  and  y-directions  are 
principal  axes  for  the  medium. 

For  an  “ofif-axis”  orientation  of  the  current  sheet,  however,  the  results  are  much  dif¬ 
ferent  and  must  be  calculated  by  using  a  simple  tensor  transformation.  Let  T**  and  Tyy 
be  the  electric  fields  produced  within  the  slab  due  to  a  current  sheet  above  the  slab.  Then 
in  the  principal  coordinate  system,  (x,y),  we  have 


so  that 


(25) 


(26) 


Let  (x',  y')  be  the  rotated  coordinate  system.  Then 

/  cosff  sintf\/x\ 
Vy'/  \-sin5  cosdMy/ 


(27) 
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Hence,  the  relationship  between  the  induced  current  and  the  exciting  current,  in  the 
(x'j  y')  coordinate  system,  is 


/  Jj/  \  /  cos  sin  /  Jj.  \ 

\Jy'  J  V  -  sin 0  COS0  J  \Jy  J 

(COS0  sind\  / aiiTxx 

—  sin  0  cos  ^  \  0 

oiiTxxCos^d -+-<Tryysin^6 
\(-anr**  +  oTyyjsin^cosd 


(28) 


(-ffiiTii  -f  aTyy)  sin  0cos0\  /  \ 

oiiTxxsin^  0  +  aTyy  cos^  0  J  \j^i^ )  ' 


The  magnitudes  of  the  currents  that  are  induced  at  the  surface  of  the  slab  are  plotted 
in  Figure  5(a),  for  a  frequency  of  10®  Hz,  and  an  =  2  x  10^,  a  =  100.  The  slab  thickness, 
zo,  is  1.27  cm  (0.5  inches),  and  the  height  of  the  current  sheet  above  the  slab,  is  2.54 
mm  (0.1  inch).  We  assume  that  the  current  sheet  carries  current  in  the  x'  direction  only; 
the  angular  variable  shown  in  Figure  5(a)  is  the  angle  between  the  x'  and  x  axes.  Now  we 
see  that,  generally,  there  is  a  component  of  current  (the  smaller  loop)  that  flows  transverse 
to  the  direction  of  the  current  sheet.  The  larger  curve  is  the  component  of  current  that 
flows  parallel  to  the  exciting  current.  For  an  isotropic  medium,  the  smaller  loop  would 
vanish,  and  the  larger  curve  would  be  an  arc  of  a  circle. 

These  results  are  in  qualitative  agreement  with  those  of  Prakash  and  Owston  [3], 
who  used  a  simple  theoretical  analysis  that  was  based  on  an  equivalent  circuit.  We  can 
extend  the  field  model  that  has  been  developed  in  this  report  to  include  slabs  that  consist 
of  graphite-epoxy  laminates  of  differing  “lay-up  order",  i.e.,  in  which  the  laminates  have 
their  principal  axes  differing  from  one-another.  In  addition,  we  can  include  the  effects  of 
non-simple  excitation  sources,  such  as  the  “horse-shoe”  eddy-current  probe  of  Prakash  and 
Owston  [3].  An  advantage  of  the  more  extensive  field-theoretic  analysis  is  that  it  more 
clearly  suggests  the  limits  of  simpler  equivalent  circuit  models. 


CIRCULAR  CURRENT  DISTRIBUTIONS 


The  circular  currents  produce  more  interesting  field  distributions  than  does  the  infinite 

^(0 

current  sheet.  Expressions  for  j  {kx,ky,z)  for  each  of  the  distributions  are  derived  in 
Appendix  C.  We  note,  now,  that  the  spectrum  of  each  current  distribution  is  no  longer 

concentrated  at  the  origin,  so  that  we  must  compute  G21  for  all  kx,ky,  as  described  in 
=  =(.i) 

Chapter  3.  When  G21  and  j  are  substituted  into  (15),  we  obtain  the  transverse  Fourier 
transform  of  the  transverse  field  vector  at  any  level,  x,  within  the  slab.  Then,  upon  taking 
the  inverse  Fourier  transform,  say  by  using  the  Fast  Fourier  Transform  (FFT)  algorithm, 
we  get  the  fields  in  physical  sp2K:e,  at  each  level,  z.  We  have  done  this  for  each  of  the  current 
distributions,  and  have  found  that  the  results  are  qualitatively  similau'  in  all  three  cases; 
hence,  we  will  display  results  of  the  filamentary  loop,  only.  The  parameters,  dimensions 
and  frequency  are  the  same  as  for  the  current  sheet;  the  radius  of  the  current  loop  is  0.5 
inch  (1.27  cm). 
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Before  going  into  the  anisotropic  problem,  we  illustrate,  in  Figure  5(b),  the  fields 
induced  into  an  isotropic  medium  (with  conductivity  2  x  10^  mhos/m)  at  a  depth  of  0.05 
inch  (0.127  cm).  The  isotropic  nature  of  the  response  is  clearly  apparent;  if  we  were  to 
look  vertically  downward  we  would  see  a  circular  response  region.  Each  pixel  is  a  square, 
whose  side  is  0.05  inch.  Thus,  the  response  region  has  a  diameter  of  about  1.0  inch,  which 
is  the  diameter  of  the  current  loop.  Therefore,  the  result  agrees  with  our  intuition. 

The  situation  in  an  anisotropic  material  is  changed  dramatically,  however.  In  this  case 
the  fibers  will  ‘guide’  the  field,  so  that  it  will  die  out  much  less  rapidly  in  the  x-direction 
(along  the  fibers)  than  in  the  y-direction.  This  is  illustrated  quite  clearly  in  Figure  5(c), 
where  the  complex  values  of  the  x-  and  y-components  of  the  electric  field  at  a  depth  of  0.05 
inch  are  shown.  The  response  region  in  this  figure  is  highly  elongated  in  the  x-direction, 
when  viewed  from  directly  above.  The  x-component  of  the  induced  electric  current  field 
is  obtained  by  multiplying  the  x-component  of  the  electric  field  by  an,  which  is  equal 
to  2  X  10^  mhos/m,  and  the  y-component  of  the  current  is  given  by  the  product  of  the 
y-component  of  electric  field  with  <722  (100  mhos/m).  Therefore,  the  eddy-currents  do  not 
flow  in  the  usual  circular  paths  of  an  unbounded  isotropic  medium,  as  suggested  by  Figure 
5(b),  but,  rather,  flow  in  highly  elongated  quasi-elliptical  paths.  The  degree  of  eccentricity 
of  the  paths  depends  upon  the  degree  of  anisotropy,  as  measured  by  the  ratio,  <711/022- 

In  many  applications  it  is  important  to  know  how  rapidly  the  induced  field  dies  out 
with  depth  into  the  slab.  In  an  anisotropic  medium  there  is  no  unique  skin-depth,  because 
the  conductivity  varies  with  direction  of  the  electric  field.  Therefore,  the  problem  must  be 
handled  numerically  in  most  cases.  We  present  in  Figure  5(d)  model  calculations  of  the 
induced  electric  field  at  a  depth  of  0.4  inches  within  the  slab.  The  excitation  source  is  the 
filamentary  current  loop  of  before,  and  the  frequency  remains  1  MHz.  Upon  comparing  this 
figure  with  Figure  5(c),  we  draw  the  following  conclusions:  the  field  magnitude  is  reduced 
by  about  an  order-of-magnitude,  the  field  is  much  more  spread  out  in  the  y-direction 
and  is  very  uniform  in  the  x-direction.  These  results  are  consistent  with  the  notion  of 
diffusion  in  isotropic  media;  they  are  an  obvious  manifestation  of  the  filtering-out  of  the 
higher  spatial  frequencies,  kx,  ky,  with  depth,  and  the  crowding  of  the  spatial-frequency 
spectrum  toward  the  origin.  In  addition,  we  note  that  the  x-component  of  the  electric  field 
dies  out  much  more  rapidly  with  depth  than  does  the  y-component.  This  is  due  to  the  fact 
that  in  the  principal  axis  coordinate  system  the  x-component  interacts  with  a  much  larger 
conductivity,  <7ii,  than  does  the  y-component  (which  interacts  with  <722-)  This  supports 
our  statement  that  there  is  no  unique  skin-depth  in  graphite-epoxy. 

Model  computations  of  the  type  presented  in  this  chapter  can  be  very  useful  in  setting 
up  eddy-current  experiments  in  graphite-epoxy,  and  interpreting  the  results. 
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CHAPTER  5 

INTEGRAL  EQUATIONS  FOR  THE  DIRECT  AND  INVERSE  PROBLEMS 
DIRECT  AND  INVERSE  PROBLEMS 

At  this  point  we  introduce  some  systems-related  ideas  that  should  make  clearer  the 
way  our  concept  of  inversion  is  used  for  nondestructive  evaluation.  Refer  to  Figure  6,  which 
shows  a  “system” ,  together  with  its  input  and  output.  In  part  (a)  of  the  figure  the  input 
is  known  auid  so  is  the  system,  and  the  output  is  to  be  determined.  This  is  the  “forward” 
or  “direct”  problem.  For  example,  the  input  could  be  current  or  voltage  source  and  the 
system,  a  coil  coupled  to  the  composite  workpiece.  The  output,  the  electromagnetic  field  or 
induced  eddy-current  within  the  workpiece,  can  be  directly  computed  in  a  straightforward 
manner  by  electromagnetic  theory;  we  have  done  this  in  Chapter  4,  by  computing  the 
eddy-current  induced  into  the  workpiece  as  a  function  of  the  orientation  of  the  exciting 
current-sheet. 

In  part  (b)  the  system  and  output  are  known,  and  the  input  is  to  be  determined. 
This  is  a  problem  of  communication  theory,  or  signal  detection.  One  assumes  a  catalog 
of  possible  input  signals  to  be  available,  whose  structure  and  characteristics  are  known  a 
priori]  from  the  known  output  one  estimates  the  input  signal  on  the  basis  of  the  maximum 
a  posteriori  probability  of  its  occurrence. 

This  is  the  basis  of  the  application  of  feature  extraction  and  artificial  intelligence  to 
nondestructive  evaluation.  It  is  an  example  of  another  forward  method,  and  appears  to  be 
sufficient  for  many  applications.  It  is,  however,  limited  by  both  the  large  volume  of  signal 
waveforms  that  must  be  catalogued  for  a  suitable  generalized  interpretation  and  by  the 
subjective  comparisons  made  by  the  interpreter.  The  method  also  gives  little  indication  of 
the  sensitivity  of  the  solution  to  possible  errors  in  the  data  and  the  degree  of  non-uniqueness 
associated  with  the  chosen  model. 

In  Figure  6(c)  both  the  input  and  output  are  known  and  the  system  is  unknown. 
The  input  could  be  a  known  probing  signal  and  the  output,  the  measured  response  to  the 
probe.  The  object  is  to  determine  the  nature  of  the'  system. 

This  is  an  example  of  system  identification,  or  pau-ameter  estimation,  where  “param¬ 
eter”  refers  to  certain  parameters  of  the  unknown  system.  In  the  sense  that  problems  (a) 
and  (b)  are  direct,  problem  (c)  is  the  “indirect”  or  “inverse”  problem,  and  is  the  problem 
attacked  in  this  research  effort.  What  one  seeks  in  this  problem  is  a  model-system  that, 
when  operating  on  the  input,  produces  a  model-output  that  is,  in  some  sense,  an  optimum 
estimate  of  the  known  output. 

In  Figure  7  we  show  a  current-sheet  that  excites  a  region  of  the  composite  workpiece, 
and  a  sensor  that  measures  the  scattered  field.  Within  the  workpiece  lies  an  anomalous 
region  (the  “fiaw”)  that  we  wish  to  reconstruct.  A  mathematical  mesh  is  defined  that 
surrounds  the  anomaly,  as  shown.  The  system,  then,  consists  of  the  current-sheet,  the 
sensor,  the  workpiece,  and  the  mesh  that  encloses  the  anomalous  region.  The  unknown 
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parameters  that  are  to  be  estimated  in  order  that  the  system  be  identified,  in  the  sense  of 
Figure  6(c)  and  the  discussion  about  it,  are  the  generalized  electrical  permittivities  that 
are  to  be  assigned  to  each  cell  of  the  mesh.  The  known  input  is  the  current  to  the  current- 
sheet,  and  the  known  outputs  are  the  field  components  at  the  sensor.  Clearly,  if  we  can 
determine  the  permittivity  map  that  is  defined  on  the  mesh,  we  will  have  reconstructed 
the  anomalous  region. 

The  method  of  solving  this  problem  is  based  on  minimizing  the  square  of  the  error 
between  the  actual  measured  data  and  that  produced  by  the  model-system,  the  model- 
output.  The  parameters  that  are  varied  to  produce  the  optimum  model,  in  the  least  squares 
sense,  are,  of  course,  the  permittivities  that  are  assigned  to  each  cell  in  the  mesh  of  Figure 
7. 

Thus,  mathematically,  we  wish  to  determine  a  set  of  unknown  parameters,  {ey},  j  = 
1, . . . ,  M ,  where  M  is  the  number  of  cells  in  the  mesh,  from  a  set  of  data,  {cj},  t  =  1, . . . ,  JV, 
where  {e,}  is  the  field  component  measured  by  the  sensor  at  a  number  of  points,  and  at  a 
number  of  different  frequencies  (if  we  are  using  a  multifrequency  approach).  The  {e^}  are 
functionally  related  to  the  {cy}  in  a  known  way: 

=/i(ci>--»fAf) 


Hence,  given  the  {cy},  we  can  calculate  the  {e^}  by  treating  this  as  a  “forward”  problem, 
in  the  sense  of  Figure  6(a). 

It  is  the  {ci},  however,  that  are  the  given  data,  so  we  must  invert  the  system,  (29), 
to  determine  the  {ey}.  We  do  this  by  minimizing  the  error  function 

. =  (30) 

«=i 

Iterative  methods  are  commonly  used  to  carry  out  the  minimization  of  (30).  The  iterative 
method  successively  improves  a  current  model,  i.e.,  a  current  estimate  of  the  {cy},  until 
the  error  measure,  (30),  is  small  and  the  parameters  are  stable  with  respect  to  reasonable 
changes  in  the  model. 

The  equations  that  define  the  forward  and  inverse  problems  are  integral  equations 
that  are  determined  by  using  electromagnetic  theory;  they  will  be  derived  next. 

INTEGRAL  RELATIONS  FOR  THE  DIRECT  AND  INVERSE  PROBLEMS 

The  detection  of  flaws,  or  anomalies,  by  means  of  eddy-currents  depends  upon  the 
fact  that  fiaws  have  different  electrical  parameters  than  the  host  material,  and,  therefore, 
the  eddy-current  flow  is  interrupted  at  the  boundary  of  the  flaw.  The  flaw,  therefore. 
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can  be  considered  to  be  an  inhomogeneity,  which  consists  of  a  tensor  permittivity,  i /  that 
is  imbedded  in  a  region  whose  tensor  permittivity,  foj  is  known  a  priori.  The  magnetic 
permeability  of  each  region  is  po.  Hence,  Maxwell’s  equations  for  the  two  regions  are: 


V  X  Eq  =  -jojpoHo 

V  X  Ho  =  jujla  ■  Eq 

known  region 

(31) (a) 

V  X  Ef  =  -jijjpoH f 

V  X  H f  =  jojtf  •  Ef 

flawed  region 

(31)(6). 

Upon  subtracting  (31)(b)  from  (3l)(a),  we  get: 

V  X  {Eq  -  Ef)  —  -joJfioiHo  - 

(32) 

V  X  {Ho  -  H f)  =  juio  •  {Eo  - 

Ef)  -h  ja;(fo  -  ^/)  •  Ef, 

where  we  have  added  and  subtracted  ju>eo  •  E /  to  get  the  final  result. 

Thus,  the  perturbation  of  the  electromagnetic  field,  Eq  —  Ef,  Hq  —  H /,  satisfies  the 
same  equation  as  the  original  electromagnetic  field  within  the  known  region,  except  for  the 
presence  of  the  term  Jo  =  Jw(lo  —  !/)•£/.  This  term,  which  is  equivalent  to  a  current 
source,  represents  the  presence  of  the  anomalous  region,  or  flaw.  It  is  important  to  note 
that  this  current  source  vanishes  off  of  the  flaw,  because  there  the  tensor  permittivities 
are  equal.  Thus,  we  say  that  the  anomalous  current  density  has  a  “compact  support”  ;  it 
occupies  a  finite  spatial  extent. 

Equation  (32)  has  the  same  form  as  (2),  with  the  anomalous  current  playing  the  role 
of  the  impressed  current.  Thus,  (32)  can  be  transformed  into  the  transverse-matrix  format, 
(11),  whose  solution  is  given  by  (15),  again  with  the  impressed  current  replaced  by  the 
anomalous  current.  The  general  form  of  the  solution  is 

(io  -  =  j  G{z\z')  ■  jjz'jdz'.  (33) 

We  suppress  the  subscript,  t,  on  the  transverse  field  variable,  e,  to  avoid  confusion  with 
other  subscripts. 

Before  going  further,  we  define: 

61,2(2)  =transverse  field  in  region  1  or  2, 
with  the  flaw  present 

60(2)  =transverse  field  with  the  flaw  absent, 
due  to  the  exciting  current  sheet. 

With  this  notation,  together  with  the  corresponding  notation  for  the  Green’s  function  in 
Chapter  3,  we  can  use  (33)  to  write  down  the  following  integral  equation  for  computing 
the  field  within  the  flawed  region,  which  is  in  region  2: 
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where 

Ja  =  -^f)  •  E2-  (34)(i>) 

The  known  field,  cq,  that  appears  in  (34)  (a)  has  already  been  calculated  in  Chapter  4. 

Equations  (34) (a), (b)  are  the  integral  equations  for  the  direct  problem:  given  the 
anomalous  permittivity,  (fo  —  f/),  compute  the  field 

The  integral  equation  for  the  inverse  problem  is  also  derived  from  (33),  except  now 
the  fields  on  the  left-hand  side  are  in  region  1,  because  they  are  measured  by  the  sensor 
(recall  Figure  7).  The  anomalous  current  is  still  in  region  2,  of  course,  so  that  we  have: 

{eo-Ei){z)=[  Gi2{z\z')  ■j^{z')dz\  (35) 

J  flaw 

where  Ja  is  defined  as  in  (34)  (b). 

Equations  (34)  and  (35)  (or,  more  precisely,  their  discretized  versions)  constitute  the 
model  system  that  was  discussed  in  (29)  and  (30).  The  rigorous  algorithm  for  using  the 
system  consists  of  first  computing  the  incident  field,  Cq,  at  the  flaw,  by  the  methods  of 
Chapter  4;  this  is  the  left-hand  side  of  (34) (a).  For  a  given  distribution  of  flaw  permittivity, 
(34)  can  be  solved  numerically  by  the  method  of  moments  [20].  The  solution  of 
(34),  62,  in  the  flawed  region  is  the  source  term  for  (35),  which  produces  the  perturbed 
field  at  the  sensor.  This  is  the  model  field  that  is  compared  with  the  measured  field  to 
determine  if  the  assumed  flaw  permittivity  is  “close”  to  the  actual  (though  unknown)  flaw 
permittivity.  The  problem  is  really  nonlinear  because  (34)  and  (35)  involve  the  product 
of  two  unknowns,  the  flaw  permittivity  and  the  electric  field  within  the  flaw.  Thus,  we 
resort  to  the  iteration  just  outlined  to  get  a  solution.  This  iteration  consists  of  a  series 
of  “direct”  problems,  which  are  the  solution  of  (34),  given  an  assumed  flaw  permittivity, 
and  “inverse”  problems,  wherein  we  determine  the  flaw  permittivity  from  (35),  given  the 
measured  fields  at  the  sensor. 

A  LINEARIZED  INVERSE  MODEL 

This  iteration  is  time  consuming,  so  we  seek  a  means  of  linearizing  the  problem.  From 
some  of  our  previous  work  in  numerical  electromagnetic  modeling  of  two-dimensional  flaws 
in  conventional  metals  [30],  together  with  laboratory  results  [31],  we  have  found  that  we 
can  approximate  62,  that  appears  in  (35),  by  the  known  field,  cq,  that  is  produced  by 
the  current  sheet,  and  that  has  already  been  computed  in  Chapter  4.  This  linearizes  the 
problem,  decouples  the  inverse  problem  from  the  direct  problem  and,  in  fact,  renders  the 
solution  of  the  direct  problem  unnecessary. 

In  the  linearized  approximation  we  have 

Ja{x,y,z)  «  ju}{Eo  -  f/)  •  Eo 
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Note  that  the  anomalous  conductivity,  a  a,  has  a  compact  support  in  the  (x,y)  plane. 

Eo  is  uniform  in  the  (x,y)  plane;  hence,  the  transverse  Fourier  transform  yields 
A:y,2)  =  &a{kx,ky,z)  •  Eq,  and  when  this  is  substituted  into  (35),  we  get  the  basic 
approximate  integral  equation  for  inversion: 


i(0  =  /  Gi2U\z')  ■  &a{z')  ■  Eo[z')dz',  (37) 

j  —Zq 

where  the  left-hand  side  is  the  perturbed  field  measured  by  the  sensor  at  2  = 

We  assume  that  the  anomalous  conductivity  tensor  has  the  same  structure  (including 
the  same  principal  ajces)  as  the  host  conductivity  tensor.  Hence,  in  the  principal  axes,  we 
have: 

[  0  0  \ 

Oa=  \  0  0  ,  (38) 

\  0  0  aa^^J 

with  d’aja  =  da33.  Our  host  material  satisfies  on  —  20,000,  <722  =  Ozz  =  100  mhos/m. 

We  know  from  our  work  with  the  current  sheet  that  an  x-directed  current  induces  an 
x-directed  electric  field,  and  a  y-directed  current  induces  a  y-directed  field  (as  long  as  we 
are  in  the  principal-axes  coordinate  system).  Hence, 

da{z')  •  Eq{z')  =  aan{z')EQi{z')ax.,  if  current  is  in  x  direction 
=  Ba-i^{z')EQ2{z')ay,  if  current  is  in  y  direction. 


DISCRETIZING  THE  INTEGRAL  EQUATION 


We  discretize  the  integral  equation  for  inversion  by  appealing  to  the  method  of  mo¬ 
ments,  together  with  a  multifrequency  approach.  In  applying  the  method  of  moments  we 
first  divide  the  slab  into  Ng  plane-parallel  regions,  and  expand  ,  and  dajj  in  terms  of 
pulse  functions  defined  on  this  partition: 


da,i  {kx,f^yjZ  )  —  dy  [kx,ky)Pj{z  ), 


where 


if  2j  ’  <  z'  <  2y 
otherwise. 


and  2y  \  Zj^^  are,  respectively,  the  bottom  and  top  of  the  jth  layer. 

When  this  is  substituted  into  (39),  and  that  result  into  (37),  we  get 
N.  +  )  ~ 

i(^)  =  ^dy*  /  Gi2(^|20  ■  axEoi{z')dz' ,  current  in  x  direction 

=  ^dy^/  Gi2{^\z')  ■  ayEo2{z')dz',  current  in  y  direction. 
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the  second.  Hence,  upon  letting  /(=  1,2, 3, 4)  stand  for  the  /th  row  of  G,  and  suppressing 

the  subscripts  on  G,  because  we  know  that  we  are  working  with  a  source  point  in  region 
2  and  field  point  in  region  1,  then  (41)  becomes 

N.  ) 

GtiiC\^')Eoi{z*)dz\  current  in  x  direction 

’ll  ^  («) 

=  Gi2{C\z')Eo2{z')dz\  current  in  y  direction. 

y=i 

This  result  is  a  discretization  of  the  unknown  conductivity  distribution.  We  proceed 
by  noting  that  the  conductivity  is  independent  of  w,  whereas  G  and  Eq  are  not.  Hence, 
we  can  get  a  system  of  equations  by  using  (42)  at  a  number,  Nf,  of  frequencies,  {u;,},  i  = 

N. 

e/(w,)  =y^ current  in  x  direction 
y=i 

(43) 

=  HljCj,  current  in  y  direction 


i  = 


The  matrices  are  defined  by: 


j  =  l,...,Nz 
/  =  1,2,3,4. 


Normally,  we  will  be  interested  in  sensing  the  magnetic  field,  and  not  the  electric. 
Thus,  /  =  3,4  in  (44),  where  ‘3’  refers  to  the  x  component,  and  ‘4’  to  the  y  component  of 
the  magnetic  field.  Thus,  we  refer  to  the  matrices  in  (44)  as  and  Hyy . 

When  the  system  matrices  of  (43)  are  nonsquare,  as  they  will  be  if  Nf  /  Ng,  then 
least-squares  methods  are  used  to  solve  the  problem.  This  model  is  called  a  “multifrequency 
model”  because  of  the  manner  of  acquiring  data.  The  frequencies  should  be  chosen  so  that 
the  system  matrix  is  reasonably  well-conditioned  in  the  mathematical  sense,  i.e.,  so  that 
the  columns  of  the  matrix  are  relatively  independent  of  each  other.  Putting  it  loosely, 
this  reduces  the  ambiguity  in  deciding  which  layer  contributes  to  the  measured  magnetic 
field.  If  we  use  a  broad  range  of  frequencies,  for  example,  the  system  matrix  will  have  a 
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structure  that  is  approximated  by  Figure  8.  The  zeroes  2U'e  not  true  zeroes  but  are  numbers 
that  are  much  smaller  than  the  x’s.  This  structure  is  due  to  the  physical  phenomenon  of 
“skin-effect” ,  wherein  high-frequency  eddy-currents  are  confined  to  the  surface  nearest  the 
exciting  source.  The  matrix  of  Figure  8,  because  it  is  roughly  triangular,  has  reasonably 
independent  columns. 

We  solve  (43)  for  each  spatial-frequency  pair,  {kx,ky),  which  gives  us  the  Fourier 
transfrom  of  the  unknown  conductivity  at  each  layer.  Then  we  take  the  inverse  Fourier 
transform  of  this  solution  vector  to  get  the  actual  distribution  of  conductivity  within  the 
anomalous  region. 


mmi 
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CHAPTER  6 

NUMERICAL  ALGORITHMS  FOR  THE  INVERSE  PROBLEM 

INTRODUCTION  TO  LEAST-SQUARES 

The  final  step  in  the  inversion  process  is  to  solve  (43) ,  which  we  rewrite  in  the  vector- 
matrix  form:  _ 

Ax  =  b,  (45) 

where  x  is  the  unknown  vector,  b  the  data  vector  and  A  the  m  x  n  system  matrix.  We 
seek  a  least-squares  solution,  as  defined  next.  Let  the  residual  vector,  r(x),  be  defined  by: 

f{x)  =  b-Ax.  (46) 

Upon  introducing  the  usual  squared-norm  notation,  we  have: 

m 

«=i 

Then  the  definition  of  a  least-squares  solution  of  (45)  is:  given  6  =  (&i, . . .  ,bm),  find  x  = 
(xi, . , .  ,Xn)  that  minimizes  l|r(2)||;  i.e.,  solve 

minj|6-A£||.  (48) 

£ 

If  X*  is  a  solution  of  (48),  then  it  is  known  that  [21] 

_ II  _ II  _ 

A  f*  =  A  (6  — Ax*)=0,  (49) 

where  the  superscript,  H,  denotes  the  complex-conjugate  transpose  (or  Hermitian  trans¬ 

pose)  of  a  matrix.  Thus, 

l”^x*  =  Z”6,  (50) 

or 

x*  =  A^6,  (51) 

=  -|-  =H  =  =H  = 

where  A  =  (A  A)  *  A  is  the  pseudoinverse  of  A. 

While  (51)  characterizes  the  mathematical  solution  of  (48),  we  don’t  actually  numeri¬ 
cally  compute  A  ,  because  of  the  potential  loss  of  precision.  For  numerical  solutions  other 
methods  are  used,  namely,  the  QR-factorization  and  singular  value  decomposition  (SVD), 
as  described  below.  These  two  methods  are  also  useful  when  A  has  less  than  full  rank,  i.e., 
when  its  columns  are  not  linearly  independent. 
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LEVENBERG-MARQUARDT  REGULARIZATION 

It  must  be  noted,  now,  that  (45)  is  a  model  equation,  whereas  the  data  vector  that 
is  actually  measured  will  contain  uncertainties  due  to  noise,  quantizing  error  in  analog-to- 
digital  conversion,  an  imprecise  model,  etc.  Due  to  this  uncertainty  in  the  data,  therefore, 
we  really  have  to  deal  with  the  following  equation: 

Ax  +  fj  =  b,  (52) 

where  fj  is  the  data-uncertainty  vector.  If  the  norm  of  ^  is  c,  then  from  (52)  we  want  the 
norm  of  the  residuals  to  satisfy: 

px  -  511  <  £.  (53) 

Since  A  is  ill-conditioned  (typical  values  of  the  condition  number  are  10^  —  10®),  the 
solution  of  (52)  will  be  unstable,  because  of  the  data  uncertainty.  To  help  smooth  the 
solution  we  impose  a  constraint  on  the  norm  of  the  solution: 

||x||  <  M.  (54) 

Inequalities  (53)  and  (54)  can  be  combined  to  yield: 

pi  -  5f  +  (jg)’||i|r  <  2£^  (55) 

Hence,  our  problem  becomes:  minimize  the  functional 

fj(i)  =  |pi-5|p+A"||lf,  (56) 

where  X  =  [(-IM)  is  the  Levenberg-Marquzwdt,  or  regularizing  parameter. 

The  vector  that  minimizes  (56)  is 


x  =  (a"a  +  A2/)-^  (57) 

Thus,  the  presence  of  the  Levenberg-Marquardt  parameter  modifies  the  pseudoinverse  of 
A  in  such  a  manner  as  to  stabilize  the  solution.  This  will  be  seen  again  when  we  discuss 
the  singular  value  decompostion  next. 

THE  SINGULAR  VALUE  DECOMPOSITION 

As  before,  we  don’t  actually  compute  a  modified  pseudoinverse;  instead,  we  work 
directly  with  A,  augmented  in  a  certain  manner.  Equation  (57)  is  the  solution  of 

___  ^  _ J| 

(A  A  +  X^J)x  =  A  5,  (58) 
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which  are  the  normal  equations  for  the  least-squares  problem  [22]: 


—  i  x  = 


where  =  denotes  equality  in  the  least-squares  sense  of  minimizing  a  norm. 

The  solution  of  (59)  is  based  on  the  singular  value  decomposition  (SVD)  [21,22]  of  A: 

A  =  t;  (  j  V  ,  (60) 

where  U{m  x  m)  and  V{n  x  n)  are  unitary,  and  S  =  diag(si,. . . ,Sn);  the  {s,}  are  the 
singular  values  of  A. 

Letting 


x  =  Vy 

=  H  _ 

g  =  U  b, 


(61)(a) 

(61)(6) 


we  get,  from  (59): 


0  P=  -- 
__  V  0 
XI 


The  XI  term  can  be  eliminated  by  using  Givens  rotations  [21,22],  with  the  result: 


0  y=  — 

I  hW 


where 


gW  ^  I  giSi/s\^\ 

^  =  diag(sS'^V..,s^^)) 


I  =  1,...  ,n 
t  =  n  -f  1, . . .  ,m 

I  =  1, . . .  ,n 


i  =  l,...,n. 


(64) (a) 

(64)  (6) 

(64) (c) 
(64) (d) 
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It  is  clear  from  (64)  (d)  that  the  Levenberg-Marquardt  parameter  accomplishes  its  stabi¬ 
lization  role  by  increzising  the  size  of  the  smallest  singular  values;  it  is  these  singular  values 
that  contribute  to  the  ill-conditioning  of  the  matrix. 

From  (63),  therefore,  we  have 


Vi  =  t  =  1, . . . ,  n  (65) 

which,  from  (61) (a),  produces  the  final  result: 

^  (66) 

»=i 

where  V  i  is  the  ith  column  vector  of  V. 

An  excellent  computer  code  for  computing  the  singular  value  decomposition  is  ZS  VDC 
in  the  LINPACK  package  [23]. 

THE  QR-FACTORIZATION 


The  advantage  of  the  singular  value  decomposition,  in  addition  to  giving  an  explicit, 
closed-form  solution,  is  that  it  produces  the  singular  values,  which  then  allows  one  to 
estimate  a  value  of  the  Levenberg-Marquardt  parameter,  by  using  (65)  and  (66),  that  will 
stabilize  the  solution.  In  addition,  one  can  then  compute  exactly  the  condition  number  of 
the  system  matrix.  A,  by  finding  the  ratio  of  the  largest  to  smallest  singular  values.  The 
disadvantage,  however,  is  that  the  computation  of  the  singular  values  and  the  matrices 
of  singular  vectors,  U,  V,  is  expensive,  in  terms  of  computer  resources  and  time.  Thus, 
we  use  another  method  to  solve  (45),  or  to  minimize  the  functional,  (56),  and  this  is  the 
QR-factorization,  which  we  now  describe,  following  [23]. 

Given  that  A  is  m  x  n  (where  m  >  n),  there  exists  a  unitary  matrix,  Q,  of  order 
m  X  m,  such  that  _ 

(67) 

where  R  is  upper  triangular.  If  we  write  Q  =  [Qi»<?2]»  where  Qi  has  n  columns,  then 

l  =  (68) 


Thus,  if  rank  (A)  =  n,  the  columns  of  Qi  form  an  orthonormal  basis  for  the  column  space 
of  A.  Now  if  A  =  [Ai,  Aj],  where  A\  has  k  columns,  and  if 


Ri2\ 
R22  J 
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where  jRn  is  k  x  k,  then 


=  H=  "11 

Q  Ai=  -- 
0 


Hence,  Q  and  iZu  give  a  QR-factorization  of  Ai.  This  truncated  decomposition  is  impor¬ 
tant  for  matrices  whose  rank  is  less  than  full,  i.e.,  for  which  rank  (A)  <  min(m,n). 

This  orthogonal  triangular ization  is  generated  by  Householder  transformations  [21,23]. 
Once  we  have  this  triangularization,  then  the  solution  to  (48)  follows: 

_ _ II  _ II  _ 

Q  f  =  Q  b-Q  Ax 


f>2  ) 

bi  ~  Rx 


Since  Q  is  unitary,  ||Q  r|l  =  |lf||.  Thus,  ||r||  is  minimized  when 


and  when  this  is  true 


fi  =  Rx, 


lr-||  =  1162 1 


Because  R  is  upper  triangular  the  solution  of  (71)  is  easily  obtained  by  back  substi¬ 
tution.  In  solving  (48)  we  apply  the  QR-factorization  to  the  matrix  A;  in  minimizing  the 
functional  (56),  we  apply  the  QR-factorization  to  the  matrix  of  the  system  (59). 

LINPACK  [23]  contains  two  useful  subroutines,  ZQRDC  and  ZQRSL,  for  solving 
least-squares  problems  by  means  of  the  QR-factorization.  ZQRDC  produces  the  QR- 
factorization  of  the  matrix  and  passes  the  result  to  ZQRSL  for  the  solution  stage. 

ANALYTIC  CONTINUATION 

We  have  already  mentioned  that  (45)  must  be  solved  for  each  spatial  frequency-pair, 
{kx,kjf),  in  order  to  take  the  inverse  Fourier  transform  of  {dy}.  The  result  is  the  conduc¬ 
tivity  distribution  within  the  flaw,  at  the  jth  layer.  Unfortunately,  the  matrix  elements  of 
A  become  ever  smaller  with  an  increase  in  (kx,ky),  thereby  rendering  A  so  ill-conditioned 
that  it  becomes  virtually  impossible  to  get  meaningful  results  for  the  lower-level  conductiv¬ 
ities  within  the  flaw.  That  is,  the  Altering  action  of  the  Levenberg-Marquardt  parameter 
forces  dj  =  0,  except  for  the  first  layer  or  two,  for  large  values  of  {kx,ky).  Or,  to  put  it 


Mm'-*  «  m  k  m  W  ft  k 
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another  way,  because  of  the  physical  process  of  field-evanescence,  which  acts  as  a  low-pass 
filter,  the  sensed  magnetic  field  at  these  high  spatial-frequencies  is  zero,  so  that,  because 
of  noise,  one  cannot  use  high  spatial-frequencies  for  reconstruction  by  solving  (45). 

High-frequency  components,  however,  are  necessary  for  high-resolution  reconstruction. 
Therefore,  it  is  necessary  to  try  to  regain  these  missing  frequencies  in  order  to  achieve  the 
desired  resolution.  That  this  attempt  at  “superresolution”  is  feasible  and  practical  (within 
limits,  of  course)  rests  upon  the  notion  of  analytic  continuation  [24,25].  We  will  state  two 
theorems  that  are  relevant  (see  [24,  p.  133]): 

Theorem  1.  The  two-dimensional  Fourier  transform  of  a  spatially  bounded  function  is 
an  analytic  function  in  the  {kx,ky)-plane. 

Theorem  2.  If  any  analytic  function  in  the  (k^,  ky)-plane  is  known  exactly  in  an  arbitrar¬ 
ily  small  (nonzero)  region  of  that  plane,  then  the  entire  function  can  be  found  (uniquely) 
by  means  of  analytic  continuation. 


1 


In  order  to  use  these  theorems  we  recall  that  a  fiaw  can  be  thought  of  as  a  space-limited 
function  g{x,y)  (or  a  function  with  compact  support).  Thus,  by  Theorem  1,  its  Fourier 
transform,  G{u,v),  is  analytic  in  the  (u,  t;)-plane.  If  we  have  only  limited  information 
about  G,  say,  only  its  values  at  low  spatial-frequencies.  Theorem  2  tells  us  that  we  can 
uniquely  extend  G  to  the  whole  (u,v)-plane.  Once  we  have  continued  G  to  the  (u,  v)-plane, 
we  can  take  an  inverse  Fourier  transform  to  recover  g. 

We  have  studied  two  approaches  to  analytic  continuation,  one  based  on  direct  matrix 
inversion  via  the  singular  value  decomposition,  and  the  other  on  an  iterative  scheme.  We 
will  discuss  only  the  latter  here. 

The  mathematical  setting  for  the  iterative  algorithm  is  projections  onto  closed  linear 
manifolds  (CLM)  in  a  Hilbert  space  [26-28].  An  important  feature  of  the  method  is  that 
a  priori  knowledge  about  g  can  be  incorporated  into  the  reconstruction  technique  in  a 
natural  way.  Such  information  includes  the  finite  support  and/or  any  constraints.  In  our 
application  we  know  that  the  conductivity  within  the  anomalous  region  is  bounded,  which 
in  our  model  implies  that  —1  <  g{x,y)  <  0,  and  we  also  have  an  estimate  of  the  support 
of  g. 

The  algorithm  that  we  used  is  listed  here: 

ALGORITHM 

Let  the  function  G{u,v)  be  given  over  a  prescribed  region  £,  and  let  J  be  the  Fourier 
transform.  Then  starting  with: 


/o(x,y)  =  7“‘[G(u,v)]; 
r  =  0; 

REPEAT 
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=  Plfr-, 

=  -Ps/,!'’; 

Pil\  =  n/i% 

F,+i=G  +  P2Fi»-, 

/,+  ,  =  p-'IPr+il; 

r  =  r  +  1; 

UNTIL  CONVERGENCE  OCCURS. 

The  important  operations  in  the  Algorithm,  in  addition  to  the  Fourier  and  inverse 
Fourier  transforms,  are  the  various  projection  operators  P,-,  which  we  now  define: 


Pi/  = 
P2F  = 


Pzf  = 


f, 

0, 

(x,y)  G  5,  S  =  support  of  /, 
otherwise. 

(73)  (a) 

G{u,u), 

(tt,t;)  €  £, 

(73)  (6) 

F{u,v), 

(u,t;)  ^  C,  where  P(u,r)  =  F[f{x,y)]. 

-1, 

if  f{z,y)  <  -1 

f{x,y), 

if-i<  /{x,y)  <0 

P3)(c) 

0, 

if  /(x,y)  >  0. 

Hence,  the  algorithm  successively  applies  the  known  properties  of  the  sought-for  solution 
to  the  initial  given  data.  The  purpose  of  the  theoretical  analyses  referred  to  in  [26-29]  is  to 
prove  that  the  algorithm  will  converge.  Numerical  experiments  suggest  various  methods 
for  speeding  the  convergence. 
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CHAPTER  7 

NUMERICAL  EXPERIMENTS 

We  have  tested  the  theory  and  numerical  algorithms,  that  were  presented  in  the 
previous  chapters  of  this  report,  with  some  numerical  experiments  that  were  performed  on 
our  computer.  The  model  analyzed  was  that  depicted  in  Figure  7,  with  the  following  data: 

thickness  of  graphite-epoxy  slab:  0.50  inches 
number  of  layers  {Ng):  10 
resolution  in  z-direction:  0.05  inches 
number  of  grid-points  in  (x,y):  64  x  64 
resolution  in  (x,y)-plane:  0.05  inches 
number  of  frequencies  (Nf):  10 

frequencies  used  (in  MHz):  1,  6,  11,  16,  21,  26,  31,  36,  45,  60 

The  values  of  the  electrical  parameters  were  those  stated  in  the  second  section  of 
Chapter  2.  The  orientation  of  the  current  is  in  the  y-direction,  which  is  transverse  to 
the  fibers.  The  x-component  of  the  magnetic  field  is  assumed  to  be  measured.  Hence, 
we  worked  with  the  matrix,  JH**',  as  defined  in  (43)  and  (44).  The  frequencies  that  were 
defined  above  worked  well  with  this  coil-sensor  orientation,  in  the  sense  that  has 
a  structure  roughly  like  that  shown  in  Figure  8,  which  is  desirable,  as  we  have  already 
discussed.  We  are  experimenting  with  other  orientations  and  frequency  ranges. 

The  flaw  consists  of  a  cylinder  of  square  cross-section  that  penetrates  vertically  the 
entire  thickness  of  the  slab;  i.e.,  the  cross-section  of  the  flaw,  in  each  of  the  ten  layers,  is 
a  square.  The  length  of  a  side  of  the  square  is  0.15  inches.  The  actual  flaw,  at  each  layer, 
looks  like  Figure  9(a). 

The  first  step  in  the  reconstruction  process  W2«  to  use  a  program  that  we  called 
RECONl.  This  program  is  our  implementation  of  the  method  described  in  the  fourth 
section  of  Chapter  6.  The  choice  of  Levenberg-Marquardt  (LM)  parameters  was  done 
empirically;  techniques  for  adaptive  selection  of  optimal  LM  parameters  are  being  studied. 
In  these  experiments,  the  parameter  values  ranged  between  10“^°  and  10“ We  point 
out  that  RECONl  computes  the  conductivities  in  the  (A;z,  ky)-plane  (Fourier  space)  for 
ail  {kx,ky).  For  large  values  of  {kx,ky),  or  high  spatial  frequencies,  the  system  equations 
become  highly  ill-conditioned  and  thus  can’t  be  expected  to  give  accurate  results  for  these 
values.  However,  for  the  small  values  of  {kx,ky),  or  low  spatial  frequencies,  the  results  are 
very  accurate.  Part  of  our  experiment  is  to  determine  how  much  low  frequency  information 
must  be  used  to  produce  fast,  accurate  reconstructions. 

In  Figure  9  we  show  the  reconstructions  of  the  flaw  at  each  layer.  These  figures  are 
labeled  with  the  name  RECON2ND  which  is  the  method  of  computing  the  inverse  FFT 
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for  all  values  of  kx  and  ky,  even  those  in  the  high  range.  As  we  can  see,  reconstruction  is 
quite  faithful  through  five  layers  and  good  from  layers  six  through  10. 

We  must  remember  that  in  order  to  use  RECON2ND,  we  must  use  RECONl  for 
all  {kx,ky).  An  alternative  is  to  use  RECONl  to  compute  the  unknowns  for  only  the 
low  frequencies,  (where  the  system  equations  are  not  highly  ill-conditioned),  and  then  use 
an  image-enhancement  technique.  The  program  LENTUY  is  our  implementation  of  the 
algorithm  described  in  the  fifth  section  of  Chapter  6. 

In  Figure  10,  we  show  the  reconstructions  generated  by  simply  taking  the  inverse  FFT 
of  the  low  frequency  data  generated  by  RECONl.  In  this  case,  we  used  — 7  <  <  7  and 
—7<ky<  7,  hence  a  ratio  of  15  to  64.  The  results  show  that  something  can  be  detected, 
but  the  reconstruction  needed  to  e  improved.  Remember,  though,  that  we  only  need  to 
produce  results  from  RECONl  for  the  low  frequencies  and,  thus,  we  save  computing  time 
and  improve  accuracy. 

Figure  11  shows  the  results  of  using  program  LENTUY  for  — 7  <  A:*  <  7,  —7<ky<  7, 
and  20  iterations.  Based  on  our  previous  experience  with  this  method,  we  believe  that  the 
number  of  iterations  can  be  reduced  to  10.  Also  for  larger  flaws,  the  values  of  kx  and  ky 
can  be  reduced.  We  see  that  this  method  improves  the  resolution  over  INVERFFT  which 
is  to  be  expected  and  that  the  reconstructions  are  faithful  through  all  10  layers  . 

Based  on  these  preliminary  experiments  and  based  on  our  previous  work  in  image 
reconstructions,  we  are  excited  by  the  results  that  have  been  generated. 
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CHAPTER  8 


COMMENTS  AND  CONCLUSIONS 

The  electrical  model  that  we  have  developed  in  this  report  is  based  on  a  continuum 
description  of  the  constitutive  relations  for  an  advanced  composite  material.  This  implies 
that  there  must  be  a  lower  limit  to  the  resolution  that  can  be  achieved  for  flaw  reconstruc¬ 
tion,  before  the  continuum  hypothesis  breaks  down  and  individual  fiber  effects  dominate. 
We  do  not  know  what  this  limit  is,  so  we  have  used  0.050  inches  as  our  resolution  in 
the  numerical  experiments.  Further  experimental  and  theoretical  analyses  are  required  to 
better  answer  this  question. 

Fiber  density  directly  affects  the  bulk  longitudinal  and  transverse  conductivities. 
Hence,  determining  (reconstructing)  these  conductivity  values  in  certain  anomalous  re¬ 
gions  allows  us  to  infer  the  status  of  the  fiber  density  within  these  regions;  it  also  informs 
us  as  to  the  possibility  of  extensive  fiber  bi-akage  within  such  regions.  Hence,  our  model 
and  algorithm  allows  us  to  draw  conclusions  about  these  classes  of  anomalies. 

In  some  cases  moisture  effects  can  be  inferred,  as  the  following  argument  suggests. 
Transverse  conduction  depends  upon  fiber-to-fiber  contact,  whereas  longitudinal  conduc¬ 
tion  does  not.  One  can  hypothesize  that  as  moisture  is  absorbed  into  the  epoxy,  the 
epoxy  swells  and  reduces  fiber  contact  [7],  thereby  reducing  the  transverse  conductivity 
but  not  the  longitudinal.  Hence,  if  we  can  reconstruct  the  transverse  conductivity  map  in 
an  anomalous  region,  as  well  as  the  longitudinal  conductivity  map,  and  if  the  transverse 
map  shows  a  significant  decrease,  whereas  the  longitudinal  map  does  not,  then  we  can  infer 
that  moisture  is  present,  and  that  there  is  no  significant  fiber  breakage,  or  fiber-density 
reduction. 

This  argument  holds  for  anisotropic  slabs,  such  ea  the  unidirectional  slab  that  we 
have  considered.  Fortunately,  our  algorithm  allows  us  to  do  such  a  reconstruction  by 
simply  orienting  the  current  sheet  along  the  fiber  direction,  and  then  transverse  to  the 
fiber  direction.  In  the  discussion  of  Chapter  5,  we  will  first  carry  out  an  inversion  using 
the  matrix  and  then  . 

The  electromagnetic  model  that  we  have  developed  is  quite  general,  in  that  it  treats 
both  isotropic  and  anisotropic  slabs.  In  addition,  it  can  be  extended  to  treat  anisotropic 
lay-ups,  i.e.,  slabs  containing  layers  arranged  in  different  directions.  Often,  however,  such 
slabs  can  be  treated  using  a  simpler  isotropic  theory,  because  the  anisotropy  effects  are 
cancelled  by  averaging  over  the  coarse  resolution  cell. 

Because  our  model  is  quite  general,  it  can  account  for  current  sources  and  induced 
currents  that  flow  normal  to  the  fiber  plane  (the  z  direction,  in  our  analysis).  Currents 
that  flow  in  the  z  direction  could  be  useful  for  distinguishing  delaminations  from  other 
anomalies,  because  delaminations  would  tend  to  perturb  these  currents  more  than  would 
other  anomalies;  more  extensive  modeling  and  experiments  during  Phase  II  will  help  clarify 
this  question,  also. 

Thus,  we  conclude  that  our  research  effort  has  met  the  specific  technical  objectives, 
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as  well  as  the  work  statement,  listed  in  Chapter  1.  Work  during  Phase  II  will  include  a 
refinement  of  the  electromagnetic  model  and  algorithm  to  account  for  more  complex  ge¬ 
ometries,  further  development  of  numerical  algorithms,  development  of  a  prototype  sensor, 
and  laboratory  experiments  with  the  prototype.  All  of  this  aims  at  the  development  of  a 
commercial  product  for  the  wide-area  inspection  and  quantification  of  advanced-composite 
structures. 
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FIGURE  1.  (a)  HOW  FIBER-TO-FIBER  CONTACT  ALLOWS  TRANSVERSE  CONDUCTION 
(REFERENCE  4). 


FIGURE  1.  (b)  A  POSSIBLE  AC  EQUIVALENT  CIRCUIT  FOR  EDDY-CURRENT  FLAW 
(REFERENCE  4). 
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FIGURE  3.  (a)  A  VECTOR  POINT-^URCE  OF  CURRENT  AT  (x',y',  z')  IN  REGION  1;  FOR 

COMPUTATION  OF  l21- 
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FIGURE  3.  (b)  A  VECTOR  POINT-OTURCE  OF  CURRENT  AT  (x',  y',  z')  IN  REGION  2;  FOR 
COMPUTATION  OF  1, 2- 


FIGURE  5.  (a)  VECTOR  CURRENT  INDUCED  WITHIN  ANISOTROPIC  SLAB,  AS  A  FUNCTION  OF 
ORIENTATION  OF  CURRENT  WITH  RESPECT  TO  FIBERS.  (THE  LARGE  CURVE 
IS  THE  MAGNITUDE  OF  THE  CURRENT  COMPONENT  PARALLEL  TO  THE  SHEET, 
AND  THE  SMALL  CURVE  IS  THE  MAGNITUDE  OF  THE  CURRENT  COMPONENT 
TRANSVERSE  TO  THE  SHEET.  FREQUENCY  =  10®  Hz;  a,,  =  2  X  10^,  a  =  100.) 
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FIGURE  5.  (b)  ELECTRIC  FIELD  INDUCED  INTO  AN  ISOTROPIC  MEDIUM,  BY  A  CIRCULAR 
FILAMENTARY  CURRENT  LOOP,  AT  A  DEPTH  OF  0.05  IN.  (FREQUENCY  = 
10®  Hz;  a,i  =  o  =  2  X  lO'*.  RE  AND  IM  DENOTE  REAL  AND  IMAGINARY 
PARTS,  RESPECTIVELY;  X  AND  Y  DENOTE  x  AND  y  COMPONENTS.) 
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FIGURE  5.  (c)  ELECTRIC  FIELD  INDUCED  INTO  AN  ANISOTROPIC  MEDIUM,  BY  A  CIRCULAR 

FILAMENTARY  CURRENT  LOOP,  AT  A  DEPTH  OF  0.05  IN.  (FREQUENCY  =  10^  Hz; 
ffll  =  2  X  lO'*,  a  =  100.  SAME  NOMENCLATURE  AS  FIGURE  5  (b).) 


mms 


NSWC  TR  85-304 


INPUT  — 
(KNOWN  I 


SYSTEM 

(KNOWN! 


OUTPUT 
(TO  BE 
COMPUTED! 


INPUT  - 

ITO  BE 
COMPUTED! 


SYSTEM 

(KNOWN! 


OUTPUT 

(KNOWN! 


INPUT  — 
(KNOWN! 


SYSTEM 
(TO  BE 
COMPUTED! 


OUTPUT 

(KNOWN! 


FIGURE  6.  SYSTEM  REPRESENTATION  OF  DIRECT  AND  INVERSE  PROBLEMS. 
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FIGURE  10.  THE  RECONSTRUCTION  OF  THE  SAME  FLAW;  USING  ONLY  THE  LOWEST 
15  SPATIAL  FREQUENCIES. 
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APPENDIX  A 

EIGENMODE  ANALYSIS  OF  ANISOTROPIC  MEDIA 


^(0 


Start  with  the  homogeneous  form  of  (11),  i.e.,  with  j  =  0.  Since  the  matrix,  5,  is 
constant  with  respect  to  z  for  a  nonstratified  medium,  we  look  for  eigenvector  solutions  of 
the  form: 

61(2)  =  eoexp(Aa),  (Al) 

where  Cq  is  a  constant  vector,  and  A  is  a  parameter  to  be  determined.  Upon  substituting 
(Al)  into  (11),  we  get  the  eigenvalue  problem: 


S  Co  =  Acq* 

The  eigenvalues  are  the  solutions  of  the  secular  equation 

det[f  -  A7]  =0. 


(A2) 


(A3) 


There  will  be  four  eigenvalues,  and  the  eigenvector  corresponding  to  each  eigenvalue  can 
be  computed  from  the  two  independent  equations  arising  out  of  the  four  simultaneous 
equations  in  (A2). 

In  the  principal  coordinate  system,  the  generalized  permittivity  tensor,  f,  is  diagonal, 
as  in  (1).  Hence,  the  only  nonzero  Kij  in  (13)  are  Kn  7^  K22  =  K33.  Upon  letting 
K22  =  Ksz  =  K,  we  get: 


kxky 


k^ 

i— 

uK' 


kl 


^14  =  -JW/io  +  J  523  =  -  j  — 


524=^^ 


uK  ’ 


S4\  =  -f  J- 


kl 


OJflO 


WHO 

k  k 

542  = 


532  =  jwK  —  j-  * 


(A4) 


WHO 


WHO 


When  these  coefficients  are  substituted  into  (A3),  and  the  determinant  expanded,  we 
obtain  a  quartic  equation  for  A: 


A^  —  cA^  "h  6  —  0, 


(A5) 


where 


o  =  —w^ho{K  +  Kii)  -|-  2ky  +  kj(l  +  KiifK) 


b  =  w^nlKKn  -  w^Ho{2Knkl  +  {K  +  /fn)*J) 
4-  klklil  +  Kn/K)  +  ktiKniK)  +  kj. 
Hence,  the  eigenvalues  (the  roots  of  (A5))  are: 

±Ai  =  ±j{w^HoKii  —  {Kii/K)kl  —  ky)^^^ 

±A3  =  ±j{w^HoK  -kl-  kiy^"^. 


(A6) 


(A7) 


A-1 
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Note  that  A3  corresponds  to  an  isotropic  medium,  and  that  the  anisotropy  is  manifest 
in  Ai.  Clearly,  if  K\i  =  K,  then  Ai  =  A3,  which  agrees  with  the  results  for  an  isotropic 
medium. 

Corresponding  to  each  eigenvalue.  A,  is  an  eigenvector  that  satisfies  (A2).  We  have 
some  liberty  in  choosing  the  two  independent  equations  that  generate  the  eigenvectors; 
hence,  there  is  some  arbitrariness  in  choosing  the  eigenvectors.  We  choose  the  following: 


‘5i4/Ai  \ 

/-5i4/Ai\ 

(  °  ^ 

f  °  ^ 

S24/A1 

0,  = 

«3  = 

1 

V4  = 

— A3/532 

1 

1  ) 

V  1  ; 

\  —SzifSzz  > 

V  —Szi/Sz2  J 

(+Ai) 

(-Ax) 

(+A3) 

(-A3) 

(Ai) 


The  second  and  fourth  vectors  are  the  two  {+)-going  modes,  and  the  first  and  third  are 
the  two  (-)-going  modes,  in  the  z-direction. 

These  results  are  used  in  computing  the  Green’s  function  in  Chapter  3. 
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APPENDIX  B 


COMPUTATION  OF  COEFFICIENTS  FOR  Gn 


In  this  appendix  we  want  to  fill  in  a  few  of  the  steps  outlined  in  Section  III  for  Gio. 
The  computations  for  G21  are  similar. 

Referring  to  Figure  3(b),  together  with  (21),  we  have  for  the  boundary  conditions  at 
z  =  0  and  z  =  —zq: 


—  0=10 
—020 
0 
1 


Uc  ' 


-he  (Bl)(a) 


g/  “2  gA.(;*'-Zo)  .  -«2  -A.(*'-*o) 

0  ^  0 


+  e'  l"  ^1'  1  +  /'  ["  ]  e-A3(«'-*o)  =  f  «20  .  ^  ^ 


0 

010 

1 

— /?20 


(Bl){6) 


The  boundary  condition  at  z  =  z'  is  the  discontinuity  induced  by  the  point  current- 
source.  For  an  x-directed  current  the  equation  is: 


7’  +« 


—  J  “2  ,  ,t  -«2  ,  -/  01 

0  +■'  0  +'  1 


\-02 


(B2)(a) 


and  for  the  y-directed  current  it  is: 


4.  ^  -“2  .  ,  01  ,f-0 

+  <i  0  +'  1  1 
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/  0  ^ 
1 

\-02J 


(52)  (6) 


Because  the  boundary  conditions  at  2  =  0  and  z  =  z'  are  independent  of  current 
direction,  it  is  best  to  work  with  them  first  and  solve  for  (c  —  /)  in  terms  of  (o,  6),  using 
(Bl)(a),  and  (c'  —  /')  in  terms  of  (s,h),  using  (Bl)(b),  The  results,  after  some  straight¬ 
forward  algebra,  are: 


e  = 


e  = 


/  = 


a(ai -aio) +feai{/32 -/92o) 

2ai  * 

^  _  o(q!i  +  aio)  +  bai{02  ~  02o) 

2ai 

-  M  -x.z' 

2/?i  * 

-a(g2<;ia-^iaao)  ^  ^ 

20t  *  ’ 

.1  _  g(Q=i  +  Qio)  +  hai{02  -  02o)  .-a, 

"  ■  2ai  " 

_  g(Q=l  -  Qio)  +  feo=l(/?2  -  ^2o)  ^X, 

2a!i 
2;8. 

'  20, 


(B3)(<.) 


e'  = 


(B3)(6) 


Now  we  substitute  these  expressions  into  the  remaining  boundary  condition  at  2  =  2', 
for  either  the  x-  or  y-directed  current  source.  For  the  x-directed  source: 


(c-c') 


+  (d-d') 


V 


and  for  the  y-directed  source: 

-ai 


(c-cO 


0 


-02)  \-02. 


.(B4)(6) 


The  solutions  of  these  equations  are: 


B-2 
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x-directed  source: 

(c-c')  =  -l/2,  (d-d')  =  -1/2,  {e-e')=0,  (/-/')=  0;  (B5)(a) 

y-directed  source: 

(c-c')  =  /32/2,  (d-d')  =  )92/2,  (e-e')  =  l/2,  (/-/')  =  1/2.  (B5)(6) 

Upon  using  (B3),  we  eliminate  (c  —  /), (c'  —  /'),  in  favor  of  {a,b,g, h): 
a(Q:i  -  aio)  +  bai{l32  -  ^o)  -  g{ai  +  aio)e^**®  -  hai{P2  -  =  vi  (B6)(a) 

a{ai  +  aio)  +  bai{P2  -  P20)  -  &(«i  -  -  hai{P2  -  -  V2  (B6)(6) 

02010 -aia20^  +  b{Pi  -  Ao)  -  -  020=10^^X3^0 

Ol  Oi  ' 

-  hipi  +  =  V3  (B6)(c) 


-a( 


0=20=10  ~  0=1020 
0=1 


)+»(A  +M  +g(— 


0=1 


-  h{Pi  -  Pio)e 


—  A3^o  _ 


V4,  (B6)(d) 


where  the  vector  on  the  right  side  is  given  by 


— oie^*** 


,  for  x-directed  source; 


(Oi/?2«^‘**  \ 

o=i/32«"^‘*’ 


,  for  y-directed  source. 


(B7)(a) 


(B7)(6) 


(B6)  and  (B7)  are  the  final  analytic  expressions;  (a,  6,  g,  h)  are  computed  numerically,  and 
the  results  used  in  the  Green’s  function.  Note  that  we  have  not  used  the  exponential  term, 
that  appears  in  (19)  and  (20).  This  is  due  to  the  fact  that  in  computing  the 
transverse  Fourier  transform  of  the  Green’s  function,  as  defined  in  (15),  we  must  divide  by 


^(0 


3 


,  which  is  the  exponential  term. 
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APPENDIX  C 

r(») 

COMPUTATION  OF  j  FOR  CIRCULAR  CURRENT  SHEETS 

If  the  sheet  lies  in  the  plane  z  =  z",  and  has  a  coil-turns  density  of  /(r),  where  r  is 
the  radial  coordinate,  then  the  current  density  is  given  by 

J^*\x,y,z)  =  Io6(z  —  z"){— Ox  sin  0  +  iiy  cos  6)  f[r),  (Cl) 

where  Iq  is  the  total  current  carried  by  the  coil.  Then 

j  {kx,ky,z)  =  °  ^  JJ  {— sin  6  +  ay  cos  6)  (C2) 


Upon  transforming  to  cylindrical  coordinates: 

\kx,ky,z)  =  ^  ^  j  »'/(r)(-a*  sintf -f  ayCOS^)e^’'(**‘^®®®+*>'®‘“®^<ir 

=  ^  ^  j  K-^xSin^-f  aycos0)e^’’(**‘=°*^+*‘'*’‘“®^d(9  dr 


The  ^-integral  can  be  easily  calculated  by  first  transforming  into  cylindrical  coordi¬ 
nates  in  Fourier-space: 

kx  =  kr  cos  <f> 

ky  =  kr  sin  <i>  (C4) 

kr  =  (kl-\-kl)'l\ 

Thus,  the  integrand  becomes  which,  according  to  a  well-known  identity 

involving  Bessel  functions,  is 

OO 

gy*,rcos(tf-^)  ^  -f  2  J^y''Jfc(fc,r)cos  k{0  -  <f>).  (C5) 

fc=i 

Only  the  first  term  survives  the  integral  over  27r  radians,  so  that  (C3)  becomes 


J  ikx,ky,z)  = 


-jIo6{z-z")  d 


+  fy(r)Mk,r)dr. 
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The  final  integral  is  the  Bessel  transform  of  the  coil-turns  density,  /(r).  This  transform 
can  be  easily  computed  for  a  number  of  interesting  practical  coil  configurations.  For 
example,  if  the  coil  consists  of  a  single  filamentary  loop  of  radius  tq,  then  /(r)  =  6[r  —  ro), 
so  that  the  Bessel  transform  is  simply  roJi(A;rro).  Hence, 


^(0  jIoro6{z- z")  ^  ,,  -  ky  .  - 

j  = - — - ^Ji(A:rro)(-a*  +  ay—). 


For  a  coil  with  a  uniform  distribution  of  turns,  whose  density  is  Nc,  and  extends 
to  a  radius  of  ro,  the  Bessel  transform  reduces  to  NcJq°  rJ\{krr)dr,  which  can  only  be 
computed  numerically.  Therefor^i,  for  such  a  coil 


^(0  jNJo6{z-z"),  _  ky  _  r 

3  = - — - {-<tx-j-  +  ay—)  rJi{krr)dr. 

ZTT  Itf  Kf  Jq 


Finally,  for  a  coil  of  radius  ro,  whose  density  of  turns  linearly  increases,  i.e.,  for  which 
/(r)  =  NcV^  the  Bessel  transform  can  be  explicitly  computed,  with  the  result  that 


zM)  3 Ncrllo6{z-  z'‘) ,  _ky  J2(A:rro)^ 
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